#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/*
 * 雅可比迭代法
 * 参数一：方程组的增广矩阵
 * 参数二：增广矩阵的行数
 * 参数三：增广矩阵的列数
 * 参数四：所要求的精度
 * 返回值：各个未知量的解所在的数值的首地址
 */
double* JacobianIterativeMethod(double Matrix[3][4],int line,int column,double precision);
/*
 * 高斯赛德尔迭代法
 * 参数一：方程组的增广矩阵
 * 参数二：增广矩阵的行数
 * 参宿三：增广矩阵的列数
 * 参数四：所要求的精度
 * 返回值：各个未知量的解所在数组的首地址
 */
double* GaussSeidelIterativeMethod(double Matrix[3][4],int line,int column,double precision);


//用于打印矩阵
void printMatrix(double Matrix[3][4],int line,int column);


int main() {

    double MatrixTest[3][4]={
            {10,-1,-2,7.2},
            {-1,10,-2,8.3},
            {-1,-1,5,4.2}
    };

    double Matrix01[3][4]={
            {50,10,1,73},
            {1,20,-1,38},
            {1,-1,20,59}
    };

    double *Ans;
//    Ans = JacobianIterativeMethod(MatrixTest,3,4,0.01);
//    for (int i = 0; i < 3; ++i) {
//        printf("Ans %d'th is %lf\n",i+1,Ans[i]);
//    }
//
//    printf("\n");
//
//    Ans = GaussSeidelIterativeMethod(MatrixTest,3,4,0.01);
//    for (int i = 0; i < 3; ++i) {
//        printf("Ans %d'th is %lf\n",i+1,Ans[i]);
//    }



    Ans = JacobianIterativeMethod(Matrix01,3,4,0.0000001);
    for (int i = 0; i < 3; ++i) {
        printf("Ans %d'th is %lf\n",i+1,Ans[i]);
    }

    printf("\n");

    Ans = GaussSeidelIterativeMethod(Matrix01,3,4,0.0000001);
    for (int i = 0; i < 3; ++i) {
        printf("Ans %d'th is %lf\n",i+1,Ans[i]);
    }

    return 0;
}

double* JacobianIterativeMethod(double Matrix[3][4],int line,int column,double precision){

    printf("Matrix is:\n"); //打印原矩阵
    printMatrix(Matrix,line,column);

    double JacobiMatrix[3][4];  //创建以雅可比方式写的矩阵
    for (int i = 0; i < line; ++i) {    //填补雅可比矩阵
        for (int j = 0; j < column; ++j) {      //对雅可比矩阵的每一行进行填补
              if(j==i)        //如果为主元素，则填补为0
                  JacobiMatrix[i][j]=0;
              else if(j==column-1)  //如果是最后一列元素，填补原矩阵最后一列元素除以原矩阵当前行主元素的值
                  JacobiMatrix[i][j]=Matrix[i][j]/Matrix[i][i];
              else      //其他为普通情况：填补原矩阵当前元素除以当前元素所在行的主元素的值乘-1的值
                  JacobiMatrix[i][j]=Matrix[i][j]/Matrix[i][i]*(-1);
        }
    }

    printf("JacobiMatrix is:\n");   //打印雅可比方法矩阵
    printMatrix(JacobiMatrix,3,4);

    double *NowAns,*PreAns;     //创建两个答案数组，NowAns用于保存当前答案，PreAns用于保存上一次的答案，用于计算精度
    NowAns = malloc(sizeof(double )*line);
    PreAns = malloc(sizeof(double )*line);
    memset(NowAns,0, sizeof(double )*line);     //数组初始化
    memset(PreAns,0, sizeof(double )*line);

    for (int i = 0; i < line; ++i)  //初始化NowAns数组
        NowAns[i]=JacobiMatrix[i][column-1];

    int flag =1;    //精度标记，为1则没有达到精度
    while (flag){

        for (int i = 0; i < line; ++i) {    //精度判断
            if(fabs(NowAns[i]-PreAns[i])>precision)     //要求每一位都达到精度
                break;
            if(i==line-1)   //每一位都达到了精度
                flag=0;
        }

        for (int i = 0; i < line; ++i) {    //保存当前答案并且清空当前答案数组
            PreAns[i]=NowAns[i];
            NowAns[i]=0;
        }
        for (int i = 0; i < line; ++i) {    //用雅可比方法计算第i+1个未知数的解
            for (int j = 0; j < line; ++j)
                NowAns[i]+=PreAns[j]*JacobiMatrix[i][j];
            NowAns[i]+=JacobiMatrix[i][column-1];
        }
    }
    return NowAns;  //返回答案
}

double* GaussSeidelIterativeMethod(double Matrix[3][4],int line,int column,double precision){

    //while中第二个for循环前的注释与雅可比方法类似
    printf("Matrix is:\n");
    printMatrix(Matrix,line,column);

    double GaussSeidelMatrix[3][4];
    for (int i = 0; i < line; ++i) {
        for (int j = 0; j < column; ++j) {
            if(j==i)
                GaussSeidelMatrix[i][j]=0;
            else if(j==column-1)
                GaussSeidelMatrix[i][j]=Matrix[i][j]/Matrix[i][i];
            else
                GaussSeidelMatrix[i][j]=Matrix[i][j]/Matrix[i][i]*(-1);
        }
    }

    printf("GaussSeidelMatrix is:\n");
    printMatrix(GaussSeidelMatrix,3,4);

    double *NowAns,*PreAns;
    NowAns = malloc(sizeof(double )*line);
    PreAns = malloc(sizeof(double )*line);
    memset(NowAns,0, sizeof(double )*line);
    memset(PreAns,0, sizeof(double )*line);

    for (int i = 0; i < line; ++i)
        NowAns[i]=GaussSeidelMatrix[i][column-1];

    int flag =1;
    while (flag){

        for (int i = 0; i < line; ++i) {
            if(fabs(NowAns[i]-PreAns[i])>precision)
                break;
            if(i==line-1)
                flag=0;
        }

        for (int i = 0; i < line; ++i) {    //用高斯赛德尔方法计算每一个未知数的解
            PreAns[i]=NowAns[i];    //保存当前未知数的解并且清空
            NowAns[i]=0;
            for (int j = 0; j < line; ++j)  //用高斯赛德尔方法求解
                NowAns[i]+=GaussSeidelMatrix[i][j]*NowAns[j];
            NowAns[i]+=GaussSeidelMatrix[i][column-1];
        }
    }
    return NowAns;
}

void printMatrix(double Matrix[3][4],int line ,int column){
    for (int i = 0; i < line; ++i) {
        for (int j = 0; j < column; ++j) {
            printf("%lf\t",Matrix[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}